Stability and decay of Bloch oscillations in presence of time-dependent nonlinearity 
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We consider Bloch oscillations of Bose-Einstein condensates in presence of a time-modulated s- 
wave scattering length. Generically, interaction leads to dephasing and decay of the wave packet. 
Based on a cyclic-time argument, we find — additionally to the linear Bloch oscillation and a rigid 
soliton solution — an infinite family of modulations that lead to a periodic time evolution of the wave 
packet. In order to quantitatively describe the dynamics of Bloch oscillations in presence of time- 
modulated interactions, we employ two complementary methods: collective-coordinates and the 
linear stability analysis of an extended wave packet. We provide instructive examples and address 
the question of robustness against external perturbations. 



I. INTRODUCTION 

Quantum dynamics of atomic Bose-Einstein conden- 
sates (BECs) in optical lattices, formed by counter- 
propagating laser beams [1-3], bears much resemblance 
to electron dynamics in solid-state crystals. For this rea- 
son, and due to the vanishingly small contribution of 
decoherence effects, BECs are an ideal test ground for 
quantum transport of matter waves in complex environ- 
ments. Perhaps one of the best examples are Bloch os- 
cillations (BOs), a phenomenon predicted by Zener [4] 
based on the band-structure framework established by 
Bloch [ ]: quantum particles in periodic potentials sub- 
jected to a constant force do not accelerate uniformly 
in real space, but oscillate instead. Because of defects 
and decoherence, BOs cannot be observed in conven- 
tional solids. With the mastery of ultracold atomic gases, 
however, BOs have been observed as a periodic motion 
of ensembles of ultracold atoms [G, 7] and BECs [1-3] in 
tilted optical lattices. 

The basic phenomenon of BOs can be well understood 
within a semi-classical framework. Let us consider a wave 
packet with a narrow momentum distribution in a lattice. 
If the wave packet is accelerated by a constant force — F 
(like gravity, for massive particles, or an electric field, 
for charged ones), then its momentum hk — —Ft will 
increase linearly. In a periodic potential with spatial 
period d, the dispersion relation of free particles is re- 
placed with a band-structure dispersion e n k with band 
index n and quasi-momcntum k. In the tight-binding 
description, which is appropriate for a very deep lat- 
tice, the lowest-band dispersion reads e k cx [1 — cos(fcd)]. 
Now, quasi-momentum and velocity are no longer pro- 
portional. Instead, a wave packet that is uniformly ac- 
celerated across the Brillouin zone has the group velocity 
v g oc dkt{k) oc sm(kd) = — sm^-Qt), oscillating with the 
Bloch frequency cjb = Fd/h. Consequently, the wave 
packet oscillates back and forth in real space. Also, re- 
lated coherent phenomena have been realized with ul- 
tracold atoms. For instance, when the external force is 



modulated harmonically in time, the interwell tunneling 
of a BEC can be suppressed [n], as predicted theoreti- 
cally in Ref. [9]. In addition, the simultaneous action of 
both constant and time-harmonic forces may lead to gi- 
ant matter-wave oscillations called super-BOs due to the 
beating of the usual BOs and the drive [10]. 

Because BOs rely on the coherent reflection of waves, 
they are very sensitive to any kind of dephasing gener- 
ated by interaction effects or lattice imperfections. Any 
deviation from perfect periodicity causes random scatter- 
ing of different /c-componcnts of the wave packet. Thus, 
its momentum distribution starts to broaden and the co- 
herent oscillations in real space are destroyed. This is 
the situation in crystalline solids, where the lattice spac- 
ing d is given by atomic distances, which are so short 
that electrons suffer from scattering events long before 
their quasi-momentum reaches the Brillouin-zone edge 
7r / 1 d. One way to overcome this problem is to artificially 
increase the lattice spacing, thus shortening the Bloch pe- 
riod, as achieved in semiconductor superlattices [11, 12]. 

Experiments with ultracold dilute gases offer very 
clean experimental conditions and open new possibili- 
ties. Atomic length scales are replaced by optical length 
scales. Atom-atom interactions — the main source of 
dephasing — are dominated by s-wave scattering. Many 
alkali species (e.g. 7 Li [13], 133 Cs [14, 15]) allow tuning 
the s-wave scattering length by means of a Feshbach res- 
onance [14, 16-18]. The s-wave scattering length can 
be tuned in a wide range, including a smooth crossover 
to negative values, i.e., attractive interaction. By sup- 
pressing the interaction entirely, one can observe very 
long-living BOs (up to 10 4 cycles in Ref. [15]). There are 
always residual experimental uncertainties, like a fluc- 
tuating s-wave scattering length, whose effect should be 
considered. Moreover, the scattering length can be delib- 
erately modulated in time, which opens a pathway to new 
effects and interesting spectroscopic applications. For ex- 
ample, a harmonic modulation in time can be used to 
probe the collective excitations of trapped BEC [19, 20]. 

In this work we present a detailed study of the stabil- 
ity and decay of BOs in tilted optical lattices for BECs 



2 



with an atom-atom contact interaction that is modulated 
harmonically in time. Throughout the article, we discuss 
all results in the BEC context. But we like to emphasize 
from the outset that our analysis is based on mathemati- 
cal properties of the nonlinear Schrodinger equation and 
thus applies to all physical systems governed by equa- 
tion (5) below. In particular, a very clean realization 
is provided by ID lattices of optical wave guides [21- 
23]. In previous work [24, 25], we have identified an infi- 
nite family of harmonic modulations g(t) that guarantee 
long-living BOs on the mean-field level. We have studied 
both the stable and unstable cases using, respectively, a 
collective-coordinates (CC) approach [26] as well as a lin- 
ear stability analysis within Floquet theory [27]. In this 
article, we re-derive these results in greater detail and 
extend them in several important aspects. 

The paper is organized as follows. In Sec. II we in- 
troduce the tight-binding approximation to the Gross- 
Pitaevskii description, suitable for BECs in deep opti- 
cal lattices. A numerical solution of the discrete Gross- 
Pitaevskii equation with time-harmonic atom-atom in- 
teraction shows the occurrence of stable and unstable 
dynamics of the BEC, depending on the frequency and 
phase of the modulation. In Sec. Ill we prove, within the 
smooth-envelope approximation, the existence of an infi- 
nite family of interactions leading to stable BOs, which is 
at odds with a quasistatic soliton stability criterion. We 
generalize the cyclic-time argument developed in Ref. [25] 
to cover all the solutions found by a different method in 
[24] , and derive the limit of validity of the wide-envelope 
ansatz. Next, Sec. IV recapitulates the CC approach 
of Ref. [24] and improves the physical interpretation in 
terms of the momentum variance. The impact of several 
relevant modulations of the interaction is discussed in de- 
tail. The CC approach captures satisfactorily the effects 
of time-dependent atom-atom interaction, as long as the 
wave-packet shape is essentially preserved. The decay of 
BOs under unstable modulations is described in Sec. V, 
where we develop a linear stability analysis of wide wave 
packets. Via perturbative Floquet theory, we study the 
growth of perturbations that ultimately destroy the wave 
packet, covering a wider class of unstable perturbations 
than presented in Ref. [25]. After discussing the respec- 
tive regimes of validity of our approaches, we summarize 
and conclude the article in Sec. VI. 



evolves according to the Gross-Pitaevskii equation 
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The homogeneous force F describes a uniform acceler- 
ation, for instance by gravity, along an axis, which we 
take to be the z axis. Even in absence of the force F, the 
initial wave packet is not the ground-state configuration, 
and the condensate tends to spread across the lattice. In 
the case of repulsive self-interaction g^D = 4irh 2 a s /m > 
(a s is the s-wave scattering length), this tendency is en- 
hanced. In the opposite case a s < 0, self-attraction coun- 
teracts dispersion and allows for soliton solutions [13]. 

Since the phenomenon of BO takes place only in the 
longitudinal z direction, we consider a lattice potential 
with strong transverse confinement, such that the trans- 
verse degrees of freedom remain frozen in their harmonic- 
oscillator ground state. In particular, we exclude from 
our study the regime of weak transverse confinement 
that results in stacks of pancake-shaped BECs. These 
are prone to transverse excitations in the presence of the 
time-dependent, sign-changing nonlinearity we consider 
below. For an extensive study of the excitation of trans- 
verse degrees of freedom, see Ref. [29]. Integrating out 
the transverse degrees of freedom then results in the one- 
dimensional Gross-Pitaevskii equation 
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with the effective interaction parameter <?id = 
mw_i_(73D/27r7i [30]. In order for Eq. (2) to be valid, the 
transverse oscillator energy Huj± must exceed all other 
relevant energies. 

If the potential Vio(z) = Vo sin 2 (tt z / d) , with lattice 
constant d, ampl itude V q > 0, and local oscillator fre- 
quency uju = 7ry/2Vo/m/d, is sufficiently deep, only the 
local harmonic oscillator ground state (or Wannier func- 
tion of the lowest band) of each lattice well is populated. 
The condensate is represented by the complex amplitudes 
\l/„ for occupying the lattice wells centered at z n = nd. A 
tight-binding equation of motion is found by integrating 
also over the z coordinate: 



im n = - J(* n+ i + * n _i) + Fdn^ n + g\$ n \H n . (3) 



II. MEAN-FIELD TIGHT-BINDING MODEL 

In typical BEC experiments, atomic gases are very di- 
lute, in the sense that the interparticle spacing exceeds 
the s-wave scattering length, which allows for a descrip- 
tion within Gross-Pitaevskii theory [28] . The condensate 
is initially created in a harmonic trap and then loaded 
into an optical lattice potential V(r) with transverse con- 
finement [15]. The condensate amplitude ^(r,t) then 



Neighboring sites are coupled by the tunneling matrix el- 
ement J « 47r- 1 /2.E r (Vb/£ r ) 3/4 exp(-2 x /Vb/S r ), where 
E x = h 2 ir 2 /(2md 2 ) is the recoil energy [30]. The tight- 
binding interaction parameter g — N*J muTy jli^hgyQ 
contains the total number of particles N because we 
choose to normalize the discrete wave function as 
Z~2 n \^n\ 2 = 1< Eq. (3) is valid only for very deep trans- 
verse and longitudinal trapping potentials, for which 
tkJ^hux > IImIU, where \\nWoo = max„ |#* 2 1 is the 
maximum local mean-field interaction energy. Under 
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these conditions the shape of the local wave functions 
does not depend much on the occupation [26, 31], and the 
tight-binding parameters J and g are also not affected. 

The tight-binding description (3) is equivalent to a 
single-band description and thus neglects Landau-Zener 
tunneling (LZT) to higher bands. Let us briefly discuss 
the conditions under which this approximation is valid. 
In the linear case (g = 0), LZT can be neglected if the 
band gap E gap = Vq/2 [30, 32] is large enough, more 
precisely, if [33] 

Fd£ r «£ g 2 ap . (4) 

In the case of constant interaction, the effective lattice 
potential is rescaled by a factor (1 + 4£ , int /i? r ) _1 [ ], 
reducing the band gap by the same factor. A typi- 
cal experimental value of Vq = AE T results in E r /J — 
8.55 » E int /J = .gl^nl 2 < 0.1 for typical parameters 
chosen below (see, e.g., Fig. 1). Thus, we find that 
the relative correction to the lattice potential E[ nt /E T 
is smaller than 1.2%, which does not change the previ- 
ous validity condition (4). Finally, in the case of mod- 
ulated interaction, e.g., g(t) = gocos(ujt + </>), one first 
has to ensure that the interaction energy remains smaller 
than the gap. For the same parameters, this results in 
Ei n t/J *C E'gap/J = 17.1, which is well fulfilled. Second, 
the modulation frequency should not become resonant 
with the gap. Below, we consider frequencies u of the 
same order as the Bloch frequency Fd/h. Thus, with 
Vo/E T of order one, the condition hu <C Sg ap is already 
included in Eq. (4) , which finally turns out to be the rel- 
evant condition also in the nonlinear cases. We conclude 
that LZT can be neglected in the situations considered 
for this work, and that the single-band description (3) is 
justified. 

Hereafter we take the lattice constant d and tunneling 
J as units of length and energy, respectively, and set 
h = l. Eq. (3) then takes the form 

M n = -tf„ +1 - # n _! + Fn* n + g(t)\V n \ 2 V n . (5) 

In the noninteracting case g = 0, the atomic cloud os- 
cillates with the Bloch frequency wb = 2%/Ts = F and 
amplitude xb ~ 2/F (in the chosen units), as we recall in 
Sec. Ill A below. A constant nonlinearity g ^ is known 
to rapidly dephase Bloch oscillations [15, 21, 26, 34]. Our 
notation g(t) emphasizes that we are interested in the ef- 
fects of an interaction that is modulated in time. In cold 
atomic gases, this can be realized by means of an external 
magnetic field close to an appropriate Feshbach resonance 
[14, 16-18]. In arrays of nonlinear optical wave-guides, 
which also allow for a tight-binding description like (5), 
time is equivalent to the propagation distance along the 
wave guides, and g(t) could be realized as a spatially 
modulated cubic nonlinearity [21, 35]. 

In order to provide a foretaste of the interesting effects 
such a time-dependent interaction can have, we solve 
Eq. (5) numerically by means of the fourth-order Runge- 
Kutta method for a wave packet with initial Gaussian 
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FIG. 1. (Color online) Condensate density |»I' n (t)| 2 as a func- 
tion of position and time, obtained by numerical integration 
of Eq. (5). Initially, the wave packet is at rest and has the 
Gaussian shape (6) with width <to = 10 chosen equal to the 
amplitude ib = 2/F, F — 0.2 of free BOs. Two different 
harmonically modulated nonlinearities (a) g(t) = go cos(Ft) 
and (b) g(t) = go sin(Ft) with go = 1 result in (a) stable and 
(b) unstable oscillations, respectively. 

shape 

vMO) = (2™ 2 )- 1 / 4 exp (-n 2 /4a 2 ) ■ (6) 

Figure 1 shows the condensate density |\I/ n (i)| 2 as a func- 
tion of position and time for two different harmonic inter- 
actions: (a) g(t) = goCOs(Ft) and (b) g(t) = gosm(Ft). 
In the first case (a), the Gaussian shape is preserved over 
time, and the interacting condensate performs long-living 
BOs with frequency ujb = F. In the second case (b) the 
initial shape is destroyed after a few cycles, as the wave 
function develops satellite peaks, and BOs are rapidly 
destroyed. 

These two cases differ solely by the relative phase be- 
tween the interaction modulation and the linear Bloch os- 
cillation, which defines a reference time starting at t = 
when the wave packet is at rest. Strikingly, in the stable 
case (a), the strongest repulsion g = +go coincides with 
the upper turning point of the wave packet, i.e. when the 
momentum is at the center of the Brillouin zone and the 
positive mass disperses the wave packet. The strongest 
attraction g = —go occurs at the lower turning point, 
i.e. with the momentum at the Brillouin zone edge and 
negative mass contracting the wave packet. Therefore, 
the observed stability clearly contradicts the simple qua- 
sistatic criterion, according to which stable BOs should 
occur when the nonlinearity compensates the lattice dis- 
persion, instead of adding to it [36] . 

The results shown in Fig. 1 prompt at least the fol- 
lowing questions, for which we will provide the answers: 

(i) Which are the periodic modulations g(t) that lead to 
stable BOs? In the upcoming Sec. Ill, we use symmetry 
considerations to identify a family of stable modulations. 

(ii) How does the interaction affect the shape of an os- 
cillating wave packet, both in stable and unstable cases? 
Section IV describes a variational approach in terms of 
collective coordinates, which provides first quantitative 
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answers, (iii) How robust are the stable cases against 
small experimental imperfections? In Sec. V, we develop 
a linear stability analysis for periodic perturbations using 
Floquet theory, which proves to be in excellent agreement 
with the numerics. 



III. CYCLIC-TIME SOLUTIONS FOR WIDE 
WAVE PACKETS 

In all of the following, we consider wide wave packets, 
which span many lattice sites, as obtained by adiabati- 
cally loading an extended BEC from a shallow trap into 
an optical lattice [15]. The first effect of the mean-field 
tight-binding equation (5) is to imprint a phase factor 
exp(— iFnt) onto each amplitude ty n . Such a phase fac- 
tor can be separated from a smooth envelope A(z,t) cen- 
tered on the moving reference point x(t), by the ansatz 



$ n (t) = e ip( - t)n A(n - x(t), t) e W 



The equation of motion obeyed by the envelope A(z, t) is 
found by Taylor-expanding the hopping terms of Eq. (5) 
as *, i± i = e ±ip (A ± d z A + d 2 A/2)e ipn+i <f>M , with third 
and higher-order derivatives of A assumed to be negli- 
gible (Sec. Ill D below discusses the validity of this as- 
sumption). The force term is taken care of by choosing 



p(t) = -Ft , 



(8) 



for the initial condition p(0) = 0. The first spatial deriva- 
tive d z A can be eliminated by setting 



x(t) = - [cos(Ft) - 1] 



(9) 



such that x(0) = 0, and <f>(t) = 2 sin(Ft) /F with 0(0) = 
without loss of generality. These choices describe the uni- 
form motion of the quasi-momentum across the Brillouin 
zone and the resulting BO in real space. During this 
oscillation, the envelope is found to obey the nonlinear 
Schrodinger equation (NLSE) 



dtA = - 



1 



lO, 



2m(t) 



d 2 A + g(t)\A\ 2 A, 



(10) 



with the oscillating inverse mass m(i)^ 1 = 2 cos(Ft). Be- 
fore analyzing the effect of a modulated interaction g(t) 
in Sec. IIIB, we first describe the usual linear BOs in the 
absence of interaction {g = 0), while paying particular 
attention to the internal breathing dynamics. 



even simpler case of constant mass by introducing the 
cyclic time 



v(t) = 



sin Ft 
F 



(11) 



Since dtf] = cos(Ft), the oscillating mass drops out 
of the resulting equation of motion id v A = —d 2 A for 
A(z, rj(t)) — A(z,t). This is the simplest free-particle 
Schrodinger equation, whose solution reads A^rf) = 
exp(— ik 2 rj)Ak(0) in momentum representation. The 
real-space solution at cyclic time rj is the initial wave 
packet A(z, 0) = A(z,0) propagated with the unitary 
evolution operator in position representation, which is 
a Gaussian. Under this evolution, an initial Gaussian 
envelope A(z, 0) = ^ z (0) such as (6) stays Gaussian, 



A(z,r]) 



2ira(rj) 



exp 



Ad{rjf 



(12) 



(7) The complex width o{ij) 2 = a 2 + irj is monotonic in the 



cyclic time rj. But expressed in the physical time rj(t) 
s'm(Ft) /F, the evolution is necessarily periodic, 



A(z,t) 



00 



2ira(t) 



exp 



Aa{tf 



(13) 



where <r(t) 2 — a 2 + ism(Ft)/F. This solution describes 
a wave packet centered at z = with variance 



a(t) 2 = / dzz 2 \A(z,t)\ 2 = <r 2 + 



sm(Ft) 2 



(14) 



The wave packet broadens only initially. At Ft = n/2, 
i.e. after the first quarter of the Bloch cycle, the mass 
changes sign and the time evolution of the width is re- 
versed. At the edge of the Brillouin zone Ft — ir, the 
wave packet recovers its original shape. Thus, the wave 
packet shows perfectly periodic breathing on top of the 
BO; instead of dispersing, it remains localized due to the 
combination of lattice and tilt. The relative amplitude 
of the breathing (14) is l/(Fa 2 ) 2 , which should be very 
small for the smooth-envelope equation (10) to be valid. 

The above discussion of the linear BO is based on the 
NLSE (10), which neglects higher derivatives, i.e., as- 
sumes that the wave packet is smooth and wide. As we 
will discuss in Sec. HID, the linear BO remains peri- 
odic beyond that assumption, even for very narrow wave 
packets. With decreasing width, the breathing increases 
and the real-space Bloch amplitude is reduced until it 
approaches zero [37, 38]. 



A. Linear Bloch oscillation with breathing 

In the linear case g — 0, Eq. (10) is the Schrodinger 
equation for a free particle with oscillating inverse mass 
m(i) -1 = 2 cos(Ft). This problem can be mapped to the 



B. Bloch-periodic interaction 

As shown already by the noninteracting solution (13), 
a Bloch-oscillating wave packet can display rich internal 
dynamics with initial broadening, provided that it recov- 
ers its initial state at the end of the Bloch period. To 
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begin with, we recall the cyclic-time argument developed 
in Ref. [25] and identify those modulations g(t) which are 
Bloch periodic and guarantee stable BOs. Later, in Sec. 
Ill C, we extend this reasoning to arbitrary (rational) fre- 
quency ratios. 

1. General case 

Motivated by the BO stability for g(t) oc cos(Ft) ob- 
served in Fig. 1, we consider the class of Bloch-periodic 
interactions that are a product of cos(Ft) and any func- 
tion P(rf) of the cyclic time rj = sin^i) /F alone: 

g(t) = co S (Ft)P( V (t)) . (15) 

Notably, this family includes the higher harmonics g(t) = 
go cos[(2n + l)Ft] and g(t) = go sm(2nFt) for all integer 
n (as well as all linear combinations thereof), because 
they can always be brought into the form (15) with the 
help of trigonometric identities. A(z,rj(t)) — A(z,t) then 
obeys the equation of motion 

id v A = -d 2 z A + P( v )\A\ 2 A , (16) 

which depends only on the cyclic time rj. And no matter 
the detailed form of its solution A{z,rj), since r/(t) is a 
periodic function of time, the solution A(z,t) must be 
periodic as well. Just as in the linear case discussed in 
Sec. Ill A above, the envelope time evolution over the first 
quarter of every Bloch period will be exactly reversed 
during the second quarter. 

The family of interactions (15) includes the cases 
g(t) = ±g cos(Ft), but not ±g sin(Ft), which is a first 
explanation of the strikingly different behavior exhibited 
in Figs. 1 (a) and (b). In the latter case, the equation 
of motion cannot be written in terms of r\ alone, so that 
the cyclic-time argument does not apply. Of course, this 
fact by itself does not necessarily imply that modulation 
(b) is unstable, but Sees. IV and V below will show that 
this is indeed the case. 



2. Special case: rigid soliton 

Let us for a moment consider the NLSE (10) from a 
quasistatic point of view, i.e. take mass m and interac- 
tion g as constant. In the usual case of positive mass, 
a linear wave packet disperses. This dispersion has to 
be compensated by an attractive interaction in order to 
obtain a stationary wave packet. For a negative mass 
(e.g. quasi-momentum close to the band edge), repulsive 
interaction is needed in order to prevent the wave packet 
from contracting. If mass and interaction have opposite 
signs, the NLSE (10) admits a stable soliton solution [39]: 

A(z,t) = -^= — - \ j- e~ iut , (17) 
K ' 1 y/M cosh (z/0 v ' 



with a characteristic width £ = — 2/(gm) > 0. 

Such a soliton configuration can be maintained during 
the BO if the interaction g is modulated such that its sign 
is always opposite to the sign of the mass [40]. However, 
the mere existence of a soliton configuration at all times 
is not sufficient for the preservation of the wave packet. 
If the equilibrium width £(t) = — 2/[m(t)g(t)] changes 
rapidly in time, the soliton cannot evolve adiabatically, 
and internal degrees of freedom are excited, which will 
finally destroy the soliton. Therefore, the simplest way to 
preserve a long-living wave packet is to have no internal 
dynamics at all. To this aim, the interaction parameter is 
modulated such that the equilibrium width £o is constant, 
i.e., g T (t) — — \g T \ cos(Ft) where \g T \ — 4/£ . Hereafter, 
this case will be referred to as rigid soliton. 

In fact, our previous discussion shows that the rigid 
soliton is but a special member of the more general fam- 
ily of stable solutions. The stable interaction of Fig. 1 (a), 
g(t) = +\go\ cos(Ft), enhances the breathing of the lin- 
ear BO, while the modulation g(t) = —\go\ cos(Ft) tends 
to suppress it, even causing antibreathing for \g \ > \g T \. 
Clearly, the —cos case fulfills the soliton stability criterion 
m{t)g{t) < for all times, whereas the +cos case does 
not. The preceding time-reversal argument, however, as- 
sures that both of them lead to undamped BOs — at least 
within the approximations underlying the NLSE (10). 
Thus, while the rigid-soliton criterion is sufficient for sta- 
bility, it is by no means necessary. 



C. Bloch-commensurate interaction 

The class of functions (15) covers all stable modula- 
tions g[t) that are higher harmonics of the Bloch fre- 
quency, i.e. with frequency uj = IF, I G N. Let us gen- 
eralize the cyclic-time argument to Bloch-commensurate 
interactions g(t) evolving at a frequency lj = IF jv (i.e. 
with period vT^jX) where v, I € N are coprime. The 
common period with cos(Ft) is then T = vT-q. 

We seek to factorize both the oscillating mass and in- 
teraction in the form 

cos(Fi) = f)M(ri) , (18) 
g(t) = r,P(r,) , (19) 

in terms of a suitable cyclic time r\ (t) and otherwise arbi- 
trary functions M and P. Then, the fj drops out of Eq. 
(10), which can be written in terms of the cyclic time rj 
only: 

id v A = ~ M(n)d 2 z A + P{rj) \A\ 2 A . (20) 

Its solution may be slightly more complicated than that 
of Eq. (16), but, again, A(z,t) — A(z,i](t)) is periodic 
because of the periodicity of r](t). 

In order to achieve the factorization (18) and (19), 
we observe that trigonometric identities permit writing 
cos Ft = ± sin vt = ± sin r M v {cos t), with a polynomial 
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constructed. If the interaction writes 




o n 



Ft = -p 



FIG. 2. (Color online) Time evolution scheme for a sta- 
ble modulation g(t) with frequency ratio l/v = 2/3. The 
common nodes of m(t) _1 = 2cos(Ft) and g(t) are shared 
by sirir, with r = r 3 i(i) = (Ft - 3vr/2)/3 [cf. Eq. (21)]. 
As required by Eq. (22), the interaction can be written as 
g(t) = 2ff sinrcosr = g sm(2F(t- 3T B /4)/3). Thus, the 
envelope function is a function of f] = cos r only and is conse- 
quently periodic in t. The right panel shows the momentum 
space density |\tfc_ p | 2 , obtained by numerical integration of 
Eq. (5) with initial condition (6) and o"o = 10. The points in 
time with the same value of n, like a and b, show the same 
distribution. Parameters are go = 5, F — 0.2. 



M„ of degree v — 1, in terms of one of the slower angles 
r defined via vt = Ft ± tt/2 mod 2-7T, or equivalently: 



(*) 



Ft 



(2j + l) , jeZ 



(21) 



Defining r\ = cost and M(rf) — (—\y(v/F)M v (rj) then 
realizes Eq. (18). And from Eq. (19), we conclude that 
all modulations of the form 



g(t) = sinrP(cosr) 



(22) 



with arbitrary P(j]) = —{F/v)P{rj), guarantee stable 
modulated Bloch oscillations. This is equivalent to two 
conditions [24]: 

• g(t) has two common zeros with cos(-Fi), namely t = 
t a , t a + T/2, the zeros of sinr [see Eq. (21)]. 

• g(t) is odd with respect to these points, g(to + 1') = 
—g{to — t'), and similar for to + T/2. 

Let us illustrate how to construct such a stable mod- 
ulation via the example of a g(t) that is harmonically 
modulated with frequency ratio uj/F — l/v = 2/3, as 
shown in Fig. 2. The commensurability v = 3 and the 
choice j — 1 determine t = T3i(t) = (Ft — 37r/2)/3, and 
from there the cyclic time r\ = cosr. According to Eq. 
(22), a harmonic modulation with the desired frequency 
ratio is g(t) = <?osin2T = 2go sinr cost. Then, the wave 
function A is a function of 77 only, as shown in the right 
panel of Fig. 2. 

Similarly, any harmonic modulation that is commen- 
surate with the Bloch frequency, cos[(l/is)Ft + 6], can be 



f I r 7r 

g(t) = g sin(7 T vj (t)) = go sin<^ - \Ft - -(2j + 1) 

(23) 

(with u,l,j e N), then the BO is stable. In Ref. [24], 
this result has been derived in a different but equivalent 
way. The class of functions (23) covers all frequencies 
commensurate with the Bloch frequency, if only the phase 
with respect to the BO is adjusted correctly (see Sec. VB 
and Fig. 5 below). 



D. Beyond the smooth-envelope approximation 

The cyclic-time argument discussed above is based on 
the NLSE (10), where higher-than-second-order spatial 
derivatives have been neglected. We still need to estimate 
the effects caused by the third-order derivative. 

In the linear case, each Fourier component can be 
treated separately, just as in Sec. Ill A. The third deriva- 
tive shows up in 



k sin Ft 



id t A k = k A |^cos Ft 

= + cos(Ft - <j> k )A k , (24) 

with tan</>fc = fc/3, which can be easily solved in terms 
of its own cyclic time rj k (t) = sm(Ft — <j> k ). The result 
is, of course, periodic, 



4 0} (t)=A k (0)e- i ^W , 



where 



<Pk(t) = 



F 



k 3 

k 2 sin Ft H (1- cos Ft) 



(25) 



(26) 



The initial Gaussian wave packet (6) with A k (0) = 
^/2/ny/ao exp(— fc 2 ^) restricts the relevant values of k 
to be of the order a^ 1 . Thus, corrections due to the third 
derivative are small for wide wave packets. The superpo- 
sition of all A k (t) leads to a reduction of the bare real- 
space Bloch amplitude (9) by (z) — [cos(F£)— 1]/(Fctq), 
which is indeed the leading order of the full result given 
in Ref. [37]. 

Interaction g(t) — g cos(Ft). What happens if inter- 
actions, as in Sec. MB, are included into Eq. (24)? The 



interaction term g(t) f dk'dqA k - 



A* 

1 A q- 



k ,A k i/2ir mixes all 



components, which consequently cannot be solved sep- 
arately. On the other hand, also the global cyclic- 
time argument from Sec. IIIB is broken by the factor 
(k/3) sm(Ft) from the third derivative. Thus, we expect 
a decay of the BO, even if the interaction satisfies the 
cyclic-time criterion (15). This decay should scale as the 
product of g and some wide-wave-packet parameter re- 
lated to o-q 1 . 

We estimate this decay analytically by considering 
the first-order correction AjL (t) — a k (t) exp[— i(p k (t)] to 
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Eq. (25), due to g(t) = g cos{Ft). 
grate the equation of motion 

dk'dq 



To this end we inte- 



«k 



-Hi 



go cos (Ft) 



4(0) (A{0) \* A 

A k-q\ A q -k>) A 



(0) 



2tt 



1(0) 



(27) 

We take advantage of the wide wave packet and expand 
systematically in the parameter a^ 1 ~ k. The inter- 
action integral evaluates to (27r)~ 1 / 4 (37ro'o)~ 1 / 2 e~ fe 
plus terms of second order in k and (Jq 1 - The time in- 
tegral of e llpk ^ cos(Ft) over the Bloch period is nonzero 
only due to the phase shift <pk ~ fc/3 between ipk(t) and 
cos (Ft). It is expressed in terms of the Bessel function 
of the first kind Ji(k 2 /F) w fc 2 /2F. Finally we get 



ak(T B ) 



~9ok 3 



exp 



(-?) 



(28) 



The relevant values of k are cut off by the exponential 
and scale as k oc (Jq 1 - No matter what the sign of go, 
the net growth of the first-order correction (28) after one 
Bloch period deforms the wave packet and destroys the 
periodic dynamics of the BO. This leaves us with the 
general scaling of the life time 



1 MT B )| ^ | g0 | 
Fife T B F^- 5 



(29) 



assuring very long life times for sufficiently wide wave 
packets. 

This result proves to be quite reliable, as we have 
checked by means of a direct integration of the tight- 
binding equation of motion (5) for several sets of param- 
eters. Numerically we defined the life time as the time 
when the momentum variance has doubled with respect 
to its initial value. This includes averaging over the con- 
tributions of many modes and dynamics that go already 
quite far away from the original perturbative perspec- 
tive. Still, we find this life time to scale as predicted by 
Eq. (29), with proportionality factor sa 0.2. 



IV. COLLECTIVE COORDINATES 

From Sec. Ill, we know which modulations g(t) should 
lead to stable BOs. But we wish to gain more quantita- 
tive information on how position and momentum distri- 
butions depend on the modulation. Moreover, we would 
like to describe the time evolution also in the unstable 
cases. Toward this aim, we employ the CC approach, as 
introduced in Ref. [26]. 



where ^ n and i^* n are canonically conjugate variables. 
Instead of describing all these amplitudes, we restrict 
the number of degrees of freedom and parametrize the 
dynamics of a smooth wave packet by its centroid x(t) = 
( n ) = Z~2n n \^n(t)\ 2 and variance w(t) = (Jn — x(t)] 2 ). 
One also needs their respective conjugate momenta p(t) 
and b(t), defined by their generating role —id p ^ n = n^ n 
and similarly for b. Thus, we employ the ansatz 



*»(*) = ±A 



^ipn-\-ib(n— x) 



(31) 



with an even envelope function A(u) that is normalized 
according to J du|_4(u)| 2 = 1 and J duu 2 \A(u)\ 2 = 1. 
The assumptions underlying the CC description (31) dif- 
fer slightly from the smooth-envelope ansatz, Eq. (7). 
Keeping a fixed wave-packet shape A(u) is of course more 
restrictive. On the other hand, the centroid x is now a 
free dynamical variable, which gives enhanced flexibility 
compared to the purely kinematic x(t) of Eq. (9). 

Inserting the CC ansatz (31) into the Hamiltonian (30), 
Taylor-expanding the discrete gradient to second order 
and performing the sum as a continuous integral, one 
finds the effective Hamiltonian 



H cc = Fx — 2 cosp 1 — 



2w 



(32) 



with the kinetic integral K = J du\A'(u)\ 2 and the inter- 
action integral / = (1/2) / dw|„4(u)| 4 . In Table I, these 
are given for a Gaussian and for a soliton-shaped wave 
packet. 

By construction, the CC variables obey the canonical 
equations of motion 



dH c , 



dx 



-F 



dH cc 
dp 



= 2 sinp 



1 - 



K + Ab 2 w 2 
2w 



dw 



w 



dH C( 
db 



_K-Aw 2 b 2 Ig(t) 
w 2 ^ 2w 3 / 2 

8wbcosp . 



(33) 
(34) 

(35) 
(36) 



In our study, the following initial conditions are consid- 
ered: x(0) = 0, p(0) = 0, w(0) = <7q, 6(0) = 0. Equation 
(33) shows that the driving term in all the equations is 
p(t) — —Ft, as already used in Eq. (8), and this inde- 
pendently of the interaction. In other words, the Bloch 
period is not affected by the atom-atom interaction. Fur- 
thermore, the dynamics of b(t) and w(t) is completely 



A. Equations of motion 

The equation of motion (5) derives as 
from the mean-field Hamiltonian 



dH/d^l 



-(*„+!*; + c.c.) +Fn|ig 



(30) 



TABLE I. Collective-coordinates parameters for Gaussian and 
soliton wave packets 

A(u) K I 



Gaussian 
soliton 



(27r)"3 e 
j^cosh ( 



7T U \ 1 



— (-V 
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defined by the autonomous Eqs. (35) and (36), whose so- 
lution then determines the centroid motion according to 
Eq. (34). 

In cold-atom experiments, time-of-flight images pro- 
vide the momentum distribution of the atomic cloud. 
Therefore, it is appropriate to study not only the average 
(quasi-) momentum p, but also the momentum variance, 
whose growth in time serves as a good indicator for the 
decay of the wave packet. Actually, the momentum vari- 
ance appears naturally already in the effective Hamil- 
tonian (32). Indeed, the kinetic energy contribution is 
nothing but the mean-field expectation — (e ij5 + e -1 ^) of 
(minus twice the real part of) the discrete translation op- 
erator. For a sufficiently narrow momentum distribution 
this contributes 



(cosp) s» cosp 1 — 



(Apf 



(37) 



Comparing with the second term on the right hand side 
of Eq. (32), we recognize (Ap) 2 = K/w + Ab 2 w. In pass- 
ing, we note that this implies that K = (Ax) 2 (Ap) 2 — 
Ab 2 (Ax) 4 is a constant of motion. Since the wave packet 
has a fixed shape, this quantity can only be the surface 
of the uncertainty ellipse, K = (Ax) 2 (Ap) 2 — (Axp) 2 , so 
that we can identify 4& 2 = (Axp) 2 / (Ax) 4 . 

In terms of the momentum variance r :— (Ap) 2 , the 
equation of motion (34) for the centroid looks slightly 
simpler: 



-2(1- -J sin (^) • 



(38) 



The momentum variance obeys the equation of motion 

r = AIg(t)b w - 1/2 . (39) 

In the linear case g = the momentum variance is 
a constant of motion, (Ap) 2 = tq = K/a 2 . 1 Then, 
Eq. (38) integrates to x(t) = (2/F) (1 - r /2) [cos(Ft) - 
1]. Compared to the lowest-order result, Eq. (9), the 
amplitude of the BO is found to be reduced by the 
finite momentum width, in agreement with the exact 
solution for Gaussian wave packets of arbitrary width, 
where the amplitude is reduced by a factor e~ r °/ 2 [37]. 
Furthermore, Eqs. (35) and (36) yield the exact solu- 
tion w(t) = cr 2 + AK sin 2 (Ft) / (Fo~ ) 2 , which agrees with 
Eq. (14). 



B. Interaction effects 

Let us now study the CC equations (33)-(36) in the 
interacting case for several different modulations of the 



1 The entire momentum distribution remains unchanged, up to the 
uniform translation p = —Ft across the Brillouin zone, as can 
be seen from the reasoning in Sec. HID: even if arbitraryly high 
derivatives are included in the first line of Eq. (24), the Fourier 
components do not mix; each of them acquires only a phase 
factor and |\l'j._p| 2 = |^4fc | 2 is stationary. 



interaction parameter g(t). In all of the following exam- 
ples, we compare the CC results to the full integration of 
the tight-binding equation. Furthermore, we analytically 
isolate the leading-order effects caused by the interaction. 
We take advantage of the wide-wave-packet condition, 
Fa 2 3> 1, and treat the interaction g perturbatively. To 
zeroth order in g and 1/Fa^, we have w(t) = cr 2 and 
b(t) = 0. Then, we compute the leading corrections via 
Eq. (35) and subsequently Eqs. (38) and (39), or Eq. (36). 



1. Constant interaction 

Constant interaction g(t) = go is known to cause 
momentum broadening and damping [26, 34]. From 
Eq. (35), we find the averaged linear increase b(t) — 
Igot/(2a^). Here and in the following, the overline de- 
notes coarse graining over one Bloch period. Via Eqs. 
(39) and (38), we find 



x(t) 





I 2 9?>t 2 








H 



cos(Ft) — 



1 - 



(40) 
(41) 



which is consistent with Refs. [26, 34]. In Fig. 3, we 
compare this perturbative result to the full CC predic- 
tion (33)-(36) and to the results from the integration 
of the discrete Gross-Pitaevskii equation (5). The ap- 
proximation (41) initially agrees nicely with the CC pre- 
diction (33)-(36). The zoomed envelope of the centroid 
motion, shown in the inset, reveals, however, that the CC 
approach initially (t < 15Tb) underestimates the damp- 
ing with respect to the full result. This is because the CC 
ansatz misses momentum broadening and energy losses 
due to degrees of freedom not included in the ansatz. At 
later times t > 20Tb, the CC approach overestimates the 
damping. At this time, the CC ansatz already begins to 
break down, because the wave packet loses its shape [Fig. 
1(b)]. 



2. Harmonic modulation 

Let us now use the CC ansatz to understand the ex- 
amples of the harmonic modulations of the interaction 
parameter presented in Fig. 1. 

a. Cosine modulation g(t) = g cos(Ft) [Fig. 1 (a)]: 
BOs are stable, in agreement with the cyclic-time argu- 
ment of Sec. III. The breathing of the linear BO is due to 
the first term in Eq. (35). The cosine modulation of g(t) 
in the second term enhances or suppresses this breathing. 
Indeed, the breathing amplitude in Eq. (14) gets multi- 
plied by a factor (l + Iaogo/2K), which in our example of 
Fig. 1 (a) evaluates to 2.82. Thus, the breathing induced 
by g Q is considerably stronger than the linear breathing. 
Nevertheless, both the integration of the CC equations 




FIG. 3. (Color online) Centroid motion of a Gaussian wave 
packet for constant interaction g(t) = go, obtained by numeri- 
cal integration of the full tight-binding equation of motion (5) 
(red dots) and CC approach (33)-(36) (solid blue). The 
dashed (green) line shows the upper turning points of the cen- 
troid within the perturbative result (41), x up — — (f — ro)/F. 
The inset shows these upper turning points on a larger scale. 
Parameters are F = 0.2, go = 1.0 and go = 10. 
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FIG. 4. (Color online) Width of the Bloch-oscillating wave 
packets of Fig. 1 (F — 0.2, ao = 10, go = 1). The respective 
curves show: Ax obtained from the tight-binding equation of 
motion (3) (red dots), the CC result Aa; = y/w (solid blue), 
the estimate Vw from Eq. (42) (dashed green). 



C. Range of validity 



of motion and the analytical result are in perfect agree- 
ment with the full tight-binding equation, as shown in 
the upper panel of Fig. 4. Thus, in this case, where the 
shape of the wave packet is preserved, the CC approach 
proves to be a very powerful tool. 

b. Sine modulation g(t) — gosin(Ft) [Fig. 1 (b)]: In 
the long run, the wave packet looses its shape, which can- 
not be accurately described by CCs. During the first few 
Bloch periods, however, the wave packet is still intact and 
we may use CCs to describe, for example, the width of the 
wave packet. The interaction g(t) enters via Eq. (35) and 
induces b(t) w K sm(Ft)/F<T$ -Ig cos(Ft)/2Fa%. Only 
the second term, proportional to go, has a non- vanishing 
time average after multiplying by cosp = cos(Ft) in the 
equation of motion (36) , which leads to 

MD^^o-^t- (42) 

This prediction is shown in the lower panel of Fig. 4, to- 
gether with the full solution of Eqs. (35) and (36) and the 
width extracted from the integration of the tight-binding 
equation of motion (5) [shown in Fig. 1 (b)]. The CC de- 
scription fares very well up to t 5Tb . At this time the 
smooth shape of the wave packet is lost, and deviations 
from the CC prediction occur without surprise. 

The change of the real-space variance (42) directly re- 
flects in the change of the momentum variance, r — r Q = 
(d w r)Aw ps 2KIg t/ {F&q), driving the decay of the BO. 
Note that for a modulation g(t) oc — sini^i, the momen- 
tum variance initially decreases. On the long run, how- 
ever, the BO still decays, because the wave packet does 
not maintain its shape, just as in the case of +sin.Ft. 



In the above examples, we have demonstrated that, 
under some limitations, the CC approach is capable of 
describing the principal degrees of freedom, position and 
variance of a Bloch-oscillating wave packet. Since the CC 
ansatz relies on the shape A(u) of the wave packet [Eq. 
(31)] to be conserved, only the stable situations can be 
well described on the long run [see Fig. 4 (a)]. In the 
unstable cases, CCs can describe only the initial dynam- 
ics, like the contraction of the wave packet shown in Fig. 
4 (b), but not the decay of its shape. For that task, we 
will pursue a different strategy in the following Section. 



V. LINEAR STABILITY ANALYSIS 

We have already seen [Figs. 1 (b) and 4 (b)] that in 
some cases perturbations on a length scale much shorter 
than the width of the wave packet occur, which ulti- 
mately destroys BOs. The periodically time-dependent 
mass and interaction, as appearing in Eq. (10), provide 
the source of energy for the growth of such perturbations, 
a phenomenon known as dynamical instability [41, 42]. 
In the following, we will employ Floquet theory to detect 
and to quantify the dynamical instability. 



A. Lyapunov exponents for excitations 

In order to quantitatively describe the growth of short- 
scale perturbations, we assume that the homogeneous 
background = ^/hq is dressed by small fluctuations 
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W„ £ C. We insert 



(43) 



into Eq. (5) and expand in powers of 5<fr n . To zeroth 
order, p(t) = —Ft and (p = — 2 cos Ft + n g(t). This 
solution describes an oscillating superfluid flow in a spa- 
tially homogeneous condensate, i.e. a 5-peak in momen- 
tum space that performs perfect BOs, and this no matter 
how strong the interaction or how it is modulated. BOs 
only get dephased by a combination of inhomogeneity 
and interaction, as we find by looking at the first-order 
equation of motion for <!>$„. Similarly to the procedure 
leading to Eq. (10), we describe 8$> x u) +z in a moving ref- 
erence frame x(t) — (2 / F)[cos(Ft) — 1]. This eliminates 
the first derivative of the Taylor expansion of S^ n ±i. 
Third and higher derivatives are assumed to be small 
and neglected. In the present approximation of a homo- 
geneous background, the equations of motion for different 
fluctuation Fourier components decouple, and their real 
and imaginary parts $fc = Sfc+idk then have the coupled 
equations of motion 



s/. = fc cos (Ft) d). , 

4 = - [k 2 cos(Ft) + 2n g(t)] s k 



(44) 
(45) 



All these components contribute to the momentum vari- 
ance 

( A P) 2 = ^E fc2 (l^| 2 + l^! 2 ) ' ( 46 ) 

k 

where N is the system size entering the discrete Fourier 
transformation. In unstable cases, some fluctuation with 
a certain fc-vector will possess the largest growth rate 
and therefore dominate the decay of BOs. We seek to 
determine the most unstable mode and its growth rate, 
or Lyapunov exponent, in the following. 

In the quasistatic picture, where the mass cos Ft 
I 2/7/ j and interaction go are considered constant, the 
system of coupled equations (44) and (45) can be solved 
by a Bogoliubov transformation [43, 44]: the elementary 
excitation 

c fc = y/ujjutf. s k + i^uj^/uJkdk , (47) 

defined in terms of single-particle dispersion uj® = 
fc 2 /(2m ) an d t ne Bogoliubov frequency 



= + 2 .9o«-o) 



(48) 



has the simple time evolution Ck{t) oc exp(— iuikt). The 
criterion for BO stability then is the following: If mo 
and go are of the same sign, then ujk is real for all fc, and 
the extended condensate is stable. If, however, mo and 
go have opposite signs, then imaginary frequencies occur 
for k < fc* = 2y /, n \g m \, indicating a modulational in- 
stability of the extended condensate. These modulations 



lead to the formation of bright solitons [39]. Consistent 
with this picture is that the critical wave number fc* , es- 
timated from the central density n = l/2£, relates to 
the soliton width: fc* = 2/£. But we have already seen 
in Sec. Ill B 2 that such a quasistatic stability criterion 
does not describe adequately the dynamical stability of 
Bloch oscillations with modulated interaction g(t). 

Let us then solve the time-dependent Eqs. (44) and 
(45) with harmonic g(t). These are linear equations with 
time-periodic coefficients, which makes them accessible 
for Floquet theory [27], provided the frequency of the ex- 
ternal modulation g(t) is commensurate with the Bloch 
frequency lob = F. Because the driving is periodic, in- 
tegrating the equations of motion over a single period 
T yields all information necessary for the time evolution 
over n € N periods: 



/ s k {t + nT) 
\d k (t + nT) 

The monodromy matrix 




sl(T) si(T) 
4(T) dl(T) 



(49) 



(50) 



contains the solution at time T starting from the two lin- 
early independent initial conditions {si(0) = l,d\(0) = 
0} and {^(O) = 0,d|(0) = 1}. From Eq. (49), it is clear 
that the eigenvalues of Mk determine the growth of 
the perturbations. With the help of Liouville's formula 
det(Mfc) = 1, one finds 



p± = (trAf fc /2) ± V(trM fe /2) 2 



1 



(51) 



The Lyapunov exponent Xk = T^ 1 In [max(|p^"|, |pZ"|)l 
then characterizes the exponential growth of the ampli- 
tudes Sk,dk ~ e Afct . In the following, we explore the 
consequences of this description in some particularly rel- 
evant cases. 



B. Robustness with respect to perturbations 

1. Stability map 

BOs are stable if the Lyapunov exponent vanishes 
for all fc [24]. This is indeed the case if g(t) fulfills 
Eq. (22), because the cyclic-time argument from Sec. Ill 
applies also in the present scope: All fluctuations Sk and 
dk evolve periodically in time, and thus do not grow ex- 
ponentially. 

Specifically, a modulation g(t) = go cos(ujt + S) at any 
frequency u) commensurate with the Bloch frequency F 
allows for periodic BOs, if the relative phase S is adjusted 
properly [Eq. (23)]. In order to assess how precise this 
phase needs to be adjusted, we can study the growth of 
fluctuations with the help of Floquet theory by comput- 
ing the life time (maxt Afc)^ 1 for a certain, small phase 
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- cos(wt) sin(wt) 



cos(arf) 



- sin(wt) 
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FIG. 5. BOs at commensurate harmonic modulation g(t) = 
go cos(ujt + 8). The stable cases (8o,u>o), according to Eq. (23), 
are marked in the 8-uj plane. Frequencies uj — IFjv have 
been accounted for up to v = 12. The size of the ellipses 
represents the life time in presence of a perturbation in phase, 
5 = 5 a + 10 -4 rad, and frequency, lo = ujq (1 + 1CT 4 ). The 
largest radii correspond to a life time of 5Tb or more, the 
smallest to 1Tb or less. Parameters are F = 0.2 and nogo = 1- 



shift of order 10~ 4 . For completeness, we also consider 
similarly small frequency detunings, although Floquet 
theory does not work because commensurability is bro- 
ken, such that we need to integrate the system of equa- 
tions (44) and (45) numerically. The life time is used as 
radius in the graphical representation in the 5-u> plane 
shown in Fig. 5. The stable points are arranged on a reg- 
ular pattern, where the most robust points are arranged 
on lines. With increasing denominator v, the robustness 
drops rapidly. 

We observe a remarkable asymmetry of the robust- 
ness between +cos Ft and —cos Ft. The stability of the 
+cos Ft modulation is much more robust than the sta- 
bility of —cos Ft. In the latter case, mass and interaction 
always have opposite sign and the quasistatic frequencies 
tOk are imaginary for k < fc* . Thus, the perturbations 
grow and decay exponentially, and their periodic return 
to the initial values relies solely on the cyclic-time argu- 
ment. Due to this exponential growth, the case —cos Ft 
is much more susceptible to perturbations than the case 
of +cos Ft, which has only real frequencies. In the fol- 
lowing, we will investigate this argument more quantita- 
tively. 



2. Perturbative Lyapunov exponents 



the Bogoliubov transformation (47), with Eq. (48) and 
2mo = 1, yields the solution Ck(rj) cx exp(— iu>kv)- It 
turns out to be handy to write the perturbed equation 
of motion for the excitation amplitude 7& = (cfc + c_k)/2 
with even parity (choosing the odd parity —i(ck — c_fc)/2 
yields the same result below): 



ij k = cos(Fi)wfe7fe + — n 5 (1) (t) (7* 



Ik, 



(52) 



So far, we have not made any approximation. Under 
the assumption that the perturbation only causes a weak 
growth of 7^ per period T, we now make the ansatz 
lk{t) = [7fc+7fc(0] ex P( — , i' UJ k T l(i)) f° r the first-order cor- 
rection Jk(t). In order to obtain the growth per total 
period T of the excitation 7^, we need to determine 



i^k Jo 



dtgM(t) ( 1 + 



0* 

2iujk s'm(Ft) I F 



Tk 



(53) 



Within the present approximation of a homogeneous 
background, the constant in the brackets of Eq. (53) con- 
tributes only via the zero-frequency component of 
causing a mere phase shift 7^ oc 17°, which may be 
dropped within the leading order. 

We now expand the perturbation in its frequency com- 
ponents 

= [9»icos(lT u0 (t))+g ul sm(lT v0 (t))] , (54) 



with Tyo(i) from Eq. (21), and v, I > and coprime. 
Within first-order perturbation theory, we may treat all 
contributions separately. In accordance with Eq. (23), 
the components g u i have vanishing contribution to the 
growth of excitations. The cosine contributions, inte- 
grated according to Eq. (53) over a fundamental period 
T = vTb, vanish for v ^ 1 and can be expressed in terms 
of Bessel functions J/ of the first kind for v = 1: 



m w fc F , 



(55) 



The relative growth ll + 7^ (T)/7ju depends on the 
complex phase 7^/(7°)* = e ia . For the Lyapunov ex- 
ponent, we are only interested in the fastest possible 
growth, i.e., we maximize with respect to a. From 
e XkT 1 + \ k T = 1 + max Q Re[7^ 1) (T)/7°] we then find 



We consider the Floquet problem (44) and (45) and 
search for analytical solutions at a given modulation 
g(t) = g W (t) + g (1) {t)- The component ff (0) (*) = 
go cos Ft does not destabilize the periodic time evolution, 
while g^ (t) is a perturbation term, which may contain 
several frequencies and phases. The unperturbed prob- 
lem is conveniently solved using the cyclic time (11), 
which eliminates the time dependence cos Ft. Then, 



Xk = 



(56) 



Within the first order considered here, the growth of the 
excitations is caused only by the components v = 1 of 
Eq. (54), i.e., integer multiples of the Bloch frequency. 
Other components v > 1 not fulfilling the cyclic-time 
condition (22) still cause growth of perturbations, but 
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only as a second-order effect in the small parameter g. 
The different contributions in Eq. (56) may add up quite 
differently, depending on the amplitudes gu- Conse- 
quently, the most unstable modes may be located at dif- 
ferent values of k. 



3. Off-phase perturbation 

The most prominent contribution to the Lyapunov ex- 
ponent (56) is the Bloch periodic perturbation sin(-Ft), 
as discussed in Ref. [25]. The Lyapunov exponent in the 
case of a Bloch periodic interaction g(t) = go cos(Ft) + 
gi sin(Fi) reads 



gi n — Ji 



2ujk 
F 



(57) 



with io k = yjk 2 (k 2 + 2g n ) [Eq. (48) with 2m = 1]. 
Here, we can connect to the phase perturbation of the 
modulation g(t) = ±go cos Ft, shown in the stability map 
of Fig. 5. We need to set g n = ±1 and gino = 10~ 4 . 
In the case go no = +1, the Bessel function in Eq. 
(57) oscillates rapidly as function of k, and the maxi- 
mum Lyapunov exponent is found close to the sixth ex- 
ternum of the Bessel function at k « 1.03 with the value 
A+ « 0.00035Tb 1 . In the case <7o n o = — L u k becomes 
imaginary for k < fc* = V2, while the analytic con- 
tinuation of ui^ 1 times the Bessel function remains real 
and regular. And indeed, the maximum Lyapunov expo- 
nent is found in the region of imaginary frequencies, at 
k « 1.05. There, the Lyapunov exponent A~ « 8.84Tg 1 
is much larger than Xj, which explains the asymmetry 
already observed in Fig. 5. Put differently, we find that 
the +cos modulation enhances the robustness against the 
sine perturbation, whereas the —cos modulation reduces 
it. 



4- Robustness against general noise 

Let us come back to the more general noise (54) and the 
prediction for the perturbation growth (56). We address 
two questions: To what extent can the robustness be 
enhanced in this case? Do the predictions hold in realistic 
systems where the wave packet is wide but finite? 

We thus confront the analytical result with a Gross- 
Pitaevskii integration [Eq. (5)], where the interaction g(f) 
is composed of an noise term of type (54) plus a deliberate 
modulation g(t) — go cos Ft. We consider different values 
of go, which in the clean case define 

• the linear BO, go = 0; 

• the breathing wave packet, go > 0; 

• the rigid soliton, g — g x < [see Eq. (17)]; 

• the antibreathing wave packet, go < gr- 
in Fig. 6, the momentum variance is shown, which signals 
the decay of the wave packet and destruction of BOs. In 
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FIG. 6. (Color online) Momentum variance (Ap) 2 of a 
wave packet of soliton shape (17) with initial spatial width 
£ = 66.16 (a small seed noise of 10 -3 mimics experimen- 
tal inhomogeneities) , Bloch-oscillating in a tilted lattice with 
F — 0.2. The interaction parameter is modulated as g(t) = 
g (i) + go cos(Ft) for a random perturbation g (1 \t) of type 
(54) with frequencies fi; = IF/u below and above the Bloch 
frequency, v = 8, I = 1,2, .. . v 2 . The coefficients g v i,g v i are 
random numbers with zero mean and standard deviation 0.4, 
except for g ulJ — 0.5 and g vv = (overwritten by —go)- The 
dash-dotted lines show the slope predicted by the Lyapunov 
exponent A — 2 maxt At [cf. Eq. (51)]. As go is varied from 
— 10 (antibreathing) over g T — —0.06 (rigid) and (linear) 
[nearly overlapping with g T ] to +10 (breathing), the predicted 
slope decreases and the observed life time increases. 



all cases, the momentum distribution starts to broaden 
at some time. Thereby, the breathing wave packet shows 
a much longer life time than the linear and the rigid 
case, whereas the antibreathing wave packet lives much 
shorter. As conjectured above, the +cos modulation in- 
deed stabilizes the BOs against noise of the interaction 
parameter g(t). 

In order to quantitatively verify the prediction (56), we 
examine the momentum density of the wave packets, as 
shown in Fig. 7 for the non-modulated and the breathing 
wave packet. We can compare the growth obtained from 
Eq. (51) and the analytical prediction (56) to the growth 
obtained by integrating the full Gross-Pitaevskii equation 
(5). Indeed, the growth of the dominantly growing mode 
(marked with a vertical line in Fig. 7) agrees with the 
largest predicted Lyapunov exponents. Beyond that, the 
full momentum-space distribution shows the finite width 
of the central wave packet and higher harmonics of the 
dominant mode, due to the nonlinearity of the Gross- 
Pitaevskii equation (5). 

For the momentum broadening (Fig. 6), the accu- 
rate description of the fastest growing mode is sufficient. 
Thus, Eq. (56) gives a satisfactory description of the BO 
decay. In conclusion, we have confirmed our prediction 
from above, that the +cos modulation of g(t) signifi- 
cantly stabilizes the BO. 
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git) = y (1) (t) git) = g (1) (t) + go cosFt 




k k 

FIG. 7. (Color online) Momentum-space portrait of the de- 
caying BOs in presence of a particular realization of the ran- 
dom perturbation g"'(t) of type (54), without (left) and with 
(right) an additional cosine-like modulation (go = 10). Up- 
per panel: Lyapunov exponents from Eq. (51) (solid black), 
and the perturbative result (56) (dashed black), which is com- 
posed of the contributions from the individual frequency com- 
ponents 0.5nouJk | Ji (2t0k/F)\ /cuk (thin gray). Lower panel: 
Stroboscopic plot of the fc-space density on a logarithmic color 
scale. Inset: Location and growth of the most unstable mode 
[prediction from full Floquet theory (51)]. Parameters are the 
same as in Fig. 6. Clearly, the harmonic modulation stabi- 
lizes the BOs against residual fluctuations of the interaction 
parameter. 



C. Range of validity 

The linear stability analysis presented in this section 
is based on the infinitely extended wave packet, and thus 
can only be expected to work for rather wide wave pack- 
ets, as in the example of Figs. 6 and 7. In other words, 
the assumption wide wave packet means that the excita- 
tions are well separated (in k space) from the main wave 
packet. 

Interactions g(t) with non-zero time average, like 
git) = go, act directly on the width degree of freedom, 
as discussed in Sec. IV B. In this case, the homogeneous 
stability analysis is not reliable, because the most im- 
portant degree of freedom is missed. Indeed, the mo- 
mentum width Ap = yfr increases linearly with time, as 
pointed out in subsection IV B, Eq. (40). Figure 8 (a) 
shows exactly this: In the case of constant nonlinearity, 
the central wave packet is spreading from the start and 
soon covers a large range in momentum space, such that 
one cannot consider excitations separated from the cen- 
tral wave packet. In this case, the dynamics is better de- 
scribed by the collective coordinates approach of Sec. IV. 
Both approaches must be considered as complementary 
to each other to describe a wide range of situations. 

On the contrary, Fig. 8 (b) shows that in the case of 




FIG. 8. (Color online) Stroboscopic plot of the momentum 
density |*l/fc| 2 in unstable cases, (a) g(t) = go- The immediate 
broadening of the wave packet in fc-space is not compatible 
with the assumptions of our linear stability analysis, whose 
starting point is an infinitely extended wave function in real 
space, (b) g(t) = gosin(Ft). The most unstable mode, which 
dominates the decay, is well separated from the central wave 
packet. In both cases, the original wave function is centered 
around k — with width oo — 100, F = 0.2 and go = 2.5. 

the sine perturbation, the perturbations remain well sep- 
arated from the central wave packet. In this case, and 
in the other cases with vanishing mean of g(t), the linear 
stability analysis proves to be very powerful. 

VI. CONCLUSIONS 

We have treated the problem of BOs with a time- 
dependent interaction in the mean-field framework of the 
one-dimensional tight-binding model. This description 
applies to a dilute Bose gas in a deep lattice potential 
with a strong transverse confinement, as well as to ar- 
rays of nonlinear optical wave-guides. Interestingly, the 
stability of BOs in presence of modulated interactions 
is sensitively conditioned on the relative phase between 
modulation and BO. 

Our analysis shows that already the linear BO has 
a breathing width, but its momentum-space distribu- 
tion is time-independent (up to the uniform translation 
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p = —Ft). For sufficiently wide wave packets, a cyclic- 
time argument allows identifying a class of interactions 
g(t) that lead to periodic dynamics, in spite of the inter- 
action. In these cases, both real-space and momentum- 
space distribution become time dependent, but return 
periodically to the initial state. In all other cases, the 
BO decays, with simultaneous momentum-space broad- 
ening. The broadening is either due to the broadening of 
the central (fc-space) wave packet, or due to the growth 
of excitations separated from the central wave packet. 

In order to quantitatively describe both the periodic 
cases and the decay, we have employed two complemen- 
tary methods, the collective- coordinates approach and 
the linear stability analysis of the extended wave packet. 
The collective-coordinates approach is valid as long as the 
shape of the wave packet is essentially conserved. It is 
capable of describing, on the one hand, the centroid and 
the breathing dynamics in the periodic cases and, on the 
other hand, the beginning of the decay in the unstable 
cases; for example, at constant interaction. The linear 
stability analysis of the extended wave packet is suitable 
for the quantitative description of the decay of wide wave 
packets, when the relevant excitations are well decoupled 
in fc-space from the original wave packet. Together, the 
two approaches provide a rather complete picture of the 
wave-packet dynamics. 

The most striking prediction due to the linear stabil- 
ity analysis is that a cosine modulation of the interac- 
tion (the one that enhances the breathing) makes the BO 
more robust with respect to certain perturbations. This 
works especially well for a fluctuating residual interaction 



with zero time average, as in Figs. 6 and 7. However, the 
modulation has little effect on the decay due to a finite 
offset of the interaction. The strategy for achieving long- 
living BOs is to tune the interaction to zero as well as 
possible and then employ the cosine modulation to min- 
imize the effect of residual fluctuations around zero. 

Conversely, we conjecture that the enhanced phase sen- 
sitivity of nonlinear BOs (similar to the effect of harmon- 
ically shaken lattices in the linear case [45]), may become 
useful to design high-precision matter-wave interferome- 
ters. 
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